## Load functions ##
source(file="/Users/ben/evo-dispersal/Quoll_model/Quoll_model_functions.R", echo=F)

###################################################################################
## Model Runs ##
load(file="/Users/ben/Documents/Papers/Current/EvoPVA/Quoll model R/Model inference/Posteriors halfpc ED Summary.RData")
post<-out
out<-c() #variable to take output
for (s in 1:100){
	alpha.s<-post["tf.alpha", 1]
	fsurv1.s<-post["tf.fsurv1", 1]
	fsurv2.s<-post["tf.fsurv2", 1]
	msurv.s<-post["tf.msurv", 1]
	beta.s<-post["tf.beta", 1]
	fec.s<-post["tf.fec", 1]
	prob.d.s<-post["tf.probd",1]
	temp<-small.mother(alpha=alpha.s, fsurv1=fsurv1.s, fsurv2=fsurv2.s, msurv=msurv.s, beta=beta.s, fec=fec.s, prob.d=prob.d.s)
	out<-rbind(out, temp)
}
colnames(out)<-c("pob4", "pob5", "pob6", "pob7", 
		"pob.sr4", "pob.sr5", "pob.sr6", "pob.sr7", "as4", "as5", "as6", "as7", "as.sr4", "as.sr5", "as.sr6", "as.sr7")
save(out, file="Posterior sims.RData")

out.sr<-out[,c(5:8, 13:16)]
out.sr[which(is.infinite(out.sr)==T)]<-NA
obs<-c(19, NA, NA, 766, 818, 470, 300, 45, NA, NA, 3228, 4820, 2590, 2890)
out<-out[,c(1:4, 9:12)]
m.sim<-apply(out, 2, median)
m.sim<-c(19, NA, NA, m.sim[1:4], 45, NA, NA, m.sim[5:8])
lower.sim<-apply(out, 2, quantile, prob=0.025)
lower.sim<-c(0, NA, NA, lower.sim[1:4], 0, NA, NA, lower.sim[5:8])
upper.sim<-apply(out, 2, quantile, prob=0.975)
upper.sim<-c(0, NA, NA, upper.sim[1:4], 0, NA, NA, upper.sim[5:8])

obs.sr<-s<-c(NA, NA, NA, 64/(18), 87/(55), 46/(2), 30/(4), NA, NA, NA, 134/(8), 167/(24), 84/1, 103/5)
m.sim.sr<-apply(out.sr, 2, median)
m.sim.sr<-c(NA, NA, NA, m.sim.sr[1:4], NA, NA, NA, m.sim.sr[5:8])
lower.sim.sr<-apply(out.sr, 2, quantile, prob=0.025)
lower.sim.sr<-c(NA, NA, NA, lower.sim.sr[1:4], NA, NA, NA, lower.sim.sr[5:8])
upper.sim.sr<-apply(out.sr, 2, quantile, prob=0.975)
upper.sim.sr<-c(NA, NA, NA, upper.sim.sr[1:4], NA, NA, NA, upper.sim.sr[5:8])

par(mfrow=c(2, 1))
ymax<-max(c(obs[1:7], upper.sim[1:7]), na.rm=T)
plot(2003:2009, obs[1:7], pch=19, ylab="Population size", xlab="",main="Pobassoo", bty="l", ylim=c(0, ymax))
points(2003:2009, m.sim[1:7], pch=19, col="grey50")
arrows(2003:2009, lower.sim[1:7], 2003:2009, upper.sim[1:7], length=0, col="grey50")
legend(2003, ymax, legend=c("Observed", "Simulated"), pch=19, col=c("black", "grey50"), bty="n")

ymax<-max(c(obs.sr[1:7], upper.sim.sr[1:7]), na.rm=T)
plot(2003:2009, obs.sr[1:7], pch=19, ylab="Sex ratio", xlab="",main="Pobassoo", bty="l", ylim=c(0, ymax))
points(2003:2009, m.sim.sr[1:7], pch=19, col="grey50")
arrows(2003:2009, lower.sim.sr[1:7], 2003:2009, upper.sim.sr[1:7], length=0, col="grey50")
legend(2003, ymax, legend=c("Observed", "Simulated"), pch=19, col=c("black", "grey50"), bty="n")

ymax<-max(c(obs[8:14], upper.sim[8:14]), na.rm=T)
plot(2003:2009, obs[8:14], pch=19, ylab="Population size", xlab="Year", main="Astell", bty="l", ylim=c(0, ymax))
points(2003:2009, m.sim[8:14], pch=19, col="grey50")
arrows(2003:2009, lower.sim[8:14], 2003:2009, upper.sim[8:14], length=0, col="grey50")
legend(2003, ymax, legend=c("Observed", "Simulated"), pch=19, col=c("black", "grey50"), bty="n")

ymax<-max(c(obs.sr[8:14], upper.sim.sr[8:14]), na.rm=T)
plot(2003:2009, obs.sr[8:14], pch=19, ylab="Sex ratio", xlab="",main="Astell", bty="l", ylim=c(0, ymax))
points(2003:2009, m.sim.sr[8:14], pch=19, col="grey50")
arrows(2003:2009, lower.sim.sr[8:14], 2003:2009, upper.sim.sr[8:14], length=0, col="grey50")
legend(2003, ymax, legend=c("Observed", "Simulated"), pch=19, col=c("black", "grey50"), bty="n")
